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Abstract — The behavior of complex aerospace systems is gov- 
erned by numerous parameters. For safety analysis it is important 
to understand how the system behaves with respect to these 
parameter values. In particular, understanding the boundaries 
between safe and unsafe regions is of major importance. In this 
paper, we describe a hierarchical Bayesian statistical modeling 
approach for the online detection and characterization of such 
boundaries. Our method for classification with active learning 
uses a particle filter-based model and a boundary-aware metric 
for best performance. From a library of candidate shapes 
incorporated with domain expert knowledge, the location and 
parameters of the boundaries are estimated using advanced 
Bayesian modeling techniques. The results of our boundary 
analysis are then provided in a form understandable by the 
domain expert. We illustrate our approach using a simulation 
model of a NASA neuro-adaptive flight control system, as well as 
a system for the detection of separation violations in the terminal 
airspace. 

I. Introduction 

Safety and controllability of an aircraft is a paramount 
requirement. Within a given flight envelope, the aircraft must, 
under all circumstances, perform safely and efficiently. Since 
the behavior of a complex aerospace system is governed 
by a multitude of parameters, it is extremely important to 
obtain knowledge on how the aircraft behaves with respect to 
these parameter values. In particular, the boundaries between 
safe and unsafe behavior are of major importance. A typical 
example is the stall speed v sta ii of an aircraft. If the current 
airspeed V as becomes smaller than v sta ii , the aircraft stalls and 
and is in severe danger of crashing. However, the stall speed 
is not a fixed number. It depends on a number of variables, 
most notably, weight, lift-over-drag, altitude, turning/climbing, 
or environmental effects like icing [1]. 

Figure 1 shows the safety boundaries for a commercial 
transport (taken from the final report of the ill-fated Air France 
flight AF447 [6]). When considering the safe and unsafe 
regions in Figure 1 with respect to aircraft speed and altitude, 
the stall speed clearly indicates boundaries, which separate safe 
areas where controlled flight is possible, from unsafe areas 
where the aircraft stalls. The non-linearity of the boundaries 
are due to physical laws and the specific aircraft design. 

For most complex systems, like aircraft, their behavior 
can only be obtained by system simulation. A typical ex- 
ample might be an adaptive flight control system, where the 
behavior of the aircraft control dynamically adapts towards 
counteracting aircraft damage or unexpected environmental 


Note : limitations shown here are valid for F-GZCP in the accident conditions (ISA+10, 
weight of 205 1. 28.7 % balance, engine and wing anti-ice ON) and MAX CLIMB thrust 



Fig. 1: Safety boundaries (grey) for aircraft over speed (Mach) 
and altitude different scenarios [6]. 

conditions. Also many computerized systems for air traffic 
control, e.g., for detection of safety violations, exhibit a highly 
complex behavior due to advanced and complex algorithms 
and their hybrid nature. Such systems often have to be treated 
as “black boxes”, i.e., there is no information available about 
their internal structure or underlying design. 

Safety boundaries for such systems can be determined 
by brute-force experiments: for each possible combination 
of input values or parameter values, a simulation is started 
and its outcome is labeled “safe” or “unsafe” accordingly. In 
practice, such an approach usually fails because a large number 
of parameters (and parameter values) prohibit the exhaustive 
exploration of the parameter space. Obviously care must be 
taken that no safety-boundary is missed, because “holes” in 
an operational envelope can have disastrous consequences and 
must therefore be recognized properly. For example, a software 
problem caused multiple computer crashes on-board a group 
of six F-22 Raptors when they crossed the 180th meridian of 
longitude (the international date line) [7]. 

A safety-boundary, once located, must be somehow char- 
acterized and described such that it can be understood by the 
domain expert. Techniques using universal function approxi- 
mators, like neural networks can be used to learn a boundary 
of any shape, but the results are hard to interpret by the domain 
expert. In practice, domain experts often have a good idea on 
how the boundary should look like. Based upon background 
knowledge and experience, they might provide information like 
’’this boundary should have a sphere-like structure”, but can 
give no indication on center and size of the sphere. For exam- 
ple, the stability boundary of an adaptive flight control system 
can be expressed as an e environment in the dimensions of the 


parameters e\, . . . , e n [2]. Then, the boundary shape could be 
expressed as JA \e 2 = e 2 for unknown A > 0 and e > 0. We 
believe that such background information substantially helps 
for locating multiple boundaries and determining their true 
shapes. 

In this paper, we describe an advanced hierarchical statisti- 
cal approach toward boundary detection and shape estimation. 
By using active learning techniques as well as efficient data 
structures and algorithms, we develop a framework that can 
detect boundaries with a low number of required simulation 
experiments. Shape and location parameters of boundaries 
are estimated using an advanced Bayesian approach, which 
is capable of providing high-quality feedback on the domain 
expert’s input in the presence of noise and system dynamic 
uncertainties. 

The remainder of this paper is structured as follows: 
Section II gives an overview of our active learning and shape 
detection approach. In Section III, we focus on boundary find- 
ing and present a metric for active learning that is aware of the 
boundary location. Section IV presents our Bayesian approach 
to modeling and analyzing the boundary shapes. Experiments 
with artificial data sets, a simulation of a NASA neuro-adaptive 
flight control, and a detection system for aircraft separation 
violations are discussed in Section V. Section VI concludes. 

II. Methodology Overview 
A. Algorithm Overview 

We propose a sequential method for the estimation of 
parameterized boundary shapes in high dimensional spaces. 
A dictionary of shape classes is provided by the domain ex- 
pert. Additional constraints on the parameters, e.g., parameter 
ranges and other prior information can be given. Typical exam- 
ples for such shape classes include (hyper-)surfaces, polygons, 
spheres, or ellipses. 

We represent our boundary problem as learning the re- 
sponse surface for the function /, where f{pc) = 1 + e if 
the experiment succeeds and f{x) = 0 + e otherwise for some 
small e > 0. In this representation a boundary is determined 
by points x with f{x) = 0.5. This representation allows us 
to formulate powerful methods to select the next data point, 
which is explained later. 

Our framework is depicted in Figure 2. Given an set of 
labeled data Do, an initial classifier is built, which provides an 
initial partitioning of the space and provides the information 
to estimate posteriors over given sets of data points. Then, 
candidate points (i.e., sets of input parameters) are iteratively 
selected by the algorithm and handed over to the computer 
experiment, which executes the system under consideration 
using the given parameters and returns a categorical result 
(success or failure). Since each run of the simulator requires 
substantial computational resources, the overall number of new 
data points should be kept as small as possible. By adding new 
data points, the classifier will be extended and improved with 
the main goal of identifying and characterizing the boundaries. 

Our algorithm is based upon the sequential classification 
and regression framework as given by DynaTree [8], [9]. It 
features dynamic regression trees and a sequential tree model. 
Particle learning for posterior simulation makes DynaTrees a 


good candidate for applications, where new data points are 
processed sequentially. In our architecture, the classifier is 
represented by a Dynatree at any given point in time. After 
adding a number of new data points, the current classifier is 
used to estimate a set of data points, which are close to the 
current prediction of the boundary. This is a subset of data 
points from a regular grid or a Latin hyper square, for which 
their entropy measure is high (classification representation) or 
the estimated response value is close to 0.5. 
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Fig. 2: Architecture overview 

This set of data points is then used to estimate the best 
parameters 0 for each of the the boundary shapes, together 
with a confidence interval for each of the parameters. 

The candidate point selection in this active learning algo- 
rithm can use as much information as is available at the current 
stage, for example, prior information given by the domain 
expert. It then selects a new point (i.e., set of input parameters), 
for which the label is obtained by running the system simulator. 
Next we present the individual steps in detail. 

III. Active learning and Expected Improvement 
A. Finding boundaries 

Each data point describing one simulation run (experiment) 
is defined as x = (Pi, . . . , P v ), where Pi are the input pa- 
rameter settings and the outcome o(x) E {success, failure}. 
Thus these data define a classification problem with Cm 2 
classes. Informally, a boundary can be found between regions, 
where all experiments yield success p(x = success) = 1 and 
those, where the experiments do not meet the success criterion 
p(x = failure) = 1. Therefore, we can define a point x to be 
on the boundary if p(x = success) = p(x = failure) = 0.5. 
Although this condition can easily be generalized to more than 
two classes, in this work, we will focus on C = 2. 

A common metric to characterize points on the bound- 
ary is based upon the entropy. The entropy entr = 
— X^cgci c c ^( x = c )l°gp( x = c ) becomes maximal at 
the boundary. In cases of more than two classes, [10] uses 
a BVSB (Best vs. Second Best) strategy. [11] defines a metric 
advantage as essentially adu(x) = | p(x = success) — p(x = 
failure)\. Then [11] considers points with minimal advantage 






to be close to the boundary. In the general case with more than 
two dimensions, [11] proposed to use the difference between 
the two most likely classes. 

In general, there are two basic methods: explicitly from 
knowledge of the classification function, or by treating the clas- 
sifier as a black box and finding the boundaries numerically. 
For some classifiers it is possible to find a simple parametric 
formula that describes the boundaries between groups, for 
example, LDA or SVM. Most classification functions can 
output the posterior probability of an observation belonging 
to a group. Much of the time we do not look at these, and just 
classify the point to the group with the highest probability. 

Points that are uncertain, i.e., have similar classification 
probabilities for two or more groups, suggest that the points 
are near the boundary between the two groups. For example, if 
a point is in Group 1 with probability 0.45, and in Group 2 with 
probability 0.55, then that point will be close to the boundary 
between the two groups. We can use this idea to find the 
boundaries. If we sample points throughout the design space 
we can then select only those uncertain points near boundaries. 
The thickness of the boundary can be controlled by changing 
the value, which determines whether two probabilities are 
similar to each other or not. Ideally, we would like this to be 
as small as possible so that our boundaries are accurate. Some 
classification functions do not generate posterior probabilities. 
In this case, we can use a k-nearest neighbors approach. 
Here we look at each point, and if all its neighbors are of 
the same class, then the point is not on the boundary and 
can be discarded. The advantage of this method is that it is 
completely general and can be applied to any classification 
function. The disadvantage is that it is slow (0(n 2 )), because 
it computes distances between all pairs of points to find the 
nearest neighbors. In general, finding of the boundaries faces 
the “curse of dimensionality”: as the dimensionality of the 
design space increases, the number of points required to make 
a perceivable boundary (for fitting or visualization purposes) 
increases substantially. This problem can be attacked in two 
ways, by increasing the number of points used to fill the design 
space (uniform grid or random sample), or by increasing the 
thickness of the boundary. 


B. Active Learning 

Computer simulation of a complex system like those dis- 
cussed above, is frequently used as a cost-effective means to 
study complex physical and engineering processes. It typically 
replaces a traditional mathematical model in cases where such 
models do not exist or cannot be solved analytically. 

Active learning, or sequential design of experiments 
(DOE), in the context of estimating response surfaces (in 
our case boundaries), is called adaptive sampling. Adaptive 
sampling starts with a relatively small space-filling input data, 
and then proceeds by fitting a model, estimating predictive 
uncertainty, and choosing future samples with the aim of min- 
imizing some measure of uncertainty, or trying to maximize 
information. We perform active learning with new data until 
the boundary is characterized with sufficient accuracy and 
confidence, and the whole space has been sufficiently explored 
to not miss any boundaries in the space. 


Consider an approach which maximizes the information 
gained about model parameters by selecting the location x, 
which has the greatest standard deviation in predicted output. 
This approach has been called ALM for Active Leaming- 
Mackay, and has been shown to approximate maximum ex- 
pected information designs [12]. An alternative algorithm is 
to select E 2 minimizing the expected reduction in the squared 
error averaged over the input space [13]. This method is 
called ALC for Active Leaming-Cohn. Rather than focusing 
on design points which have large predictive variance, ALC 
selects configurations that would lead to a global reduction in 
predictive variance. 

The ALM/ ALC algorithms are suitable for classification 
but not primarily for boundary detection [14]. These heuristics 
are in general not suited for the boundary-finding task because 
they do not take the specifics of the boundaries into account 
and they tend to also explore sparsely populated regions far 
away from current boundaries. 

C. Boundary Expected Improvement 

Finding a boundary between two classes can be considered 
as finding a contour with a = 0.5 in the response surface of 
the system response. Inspired by [16] and work on contour 
finding algorithms, we loosely follow [15], and define our 
heuristics by using an improvement function. In order to use 
the available resources as efficiently as possible for our con- 
tour/boundary finding task, one would ideally select candidate 
points which lie directly on the boundary, but that is unknown. 
Therefore, new trial points x are selected, which belong to 
an e-environment around the current estimated boundary. This 
means that 0.5 — e < y(x) < 0.5 + e. New data points should 
maximize the information in the vicinity of the boundary. 
Following [16] and [15], we define an improvement function 
for x as 

I(x) = e 2 (x) — min {(y(x) — 0.5) 2 , e 2 (x)} 

Here, y{x) ~ N(y(x), cr 2 (x)), and e(x) = acr(x) for a 

constant a > 0. This term defines an e-neighborhood around 
the boundary as a function of a(x). This formulation makes 
it possible to have a zero-width neighborhood around existing 
data points. For boundary sample points, I(X) will be large 
when the predicted a(x) is largest. 

The expected improvement E[I(x)\ can be calculated eas- 
ily following [15] as 

0.5+acr(:r) 

E[I(X)} = - I (y-y( X ))^(^-M^dy 

0.5—acr(x) 

+2 (y(x) - 0.5)cr 2 (a;) [<j>(z + (x)) - <j>(z-(x))\ 
+(aV0r) - (y(x) - 0.5) 2 ) [$(*+(*)) ~ *(*-(*))] , 

where z± (x) = ± a, and <j> and $ are the standard 

normal density and cumulative distribution, respectively. Each 
of these three terms are instrumental in different areas of the 
space. The first term summarizes information from regions of 
high variability within the e-band. The integration is performed 
over the e-band as e(x) = aa(x). The second term is 
concerned with areas of high variance farther away from the 
estimated boundary. Finally, the third term is active close to 



the estimated boundary. After the expected improvement has 
been calculated, the candidate point is selected as the point, 
which maximizes the expected improvement. 

IV. Boundary shape 

1 ) Notation : Suppose there are m shape classes 

with m > 1, which are parameterized by 
©i, . . . , 0 m . The task is to fit l shapes Si, . . . , S/, l > 1, 
where Si = (ii, ©i Si = and ij denotes the 

shape class for the j th shape with ij G M = {Mi, . . . , M m }. 
Several of the ij can be the same to accommodate more than 
one shape belonging to the same class. The 0^ should be 
different since we do not want to represent the same boundary 
shape twice. We also seek to determine the correct number of 
shapes l that represents the input point set X n . 

For example, we may consider the m = 2 shape classes 
Mi = hyperplane and M 2 = sphere in R d . Hyperplanes 
are represented as a\X\ + • • • + a^x^ + a^+i = 0 with 
parameter vector 0i = (ai, . . . , a^, a^+i) G R d+1 . In the 
same d-dimensional space, a sphere of radius r with center 

c = (ci, . . . , Cd) is described by {pc\ — ci) 2 H h(^d — cf) 2 = 

r 2 with parameter vector 0 2 = (c, r) G i? d+1 . 

2) A a Good Shape Set S for an Input Point Set 
X n ?: There are three conditions that specify when a shape 
set S provides a good fit to the data X n \ 

(i) Summary: each point on a shape S G S is close to 
some classifier boundary point in X n , 

(ii) Completeness : each classifier boundary point in X n is 
close to some shape point on some shape S G <S, and 

(iii) Minimality : the shapes in S are as different from one 
another as possible. 

Condition (i) encourages each shape S G S to be a good 
summary of one of the parts of the boundary of classifier P n . 
That is, the points of a shape should lie along high entropy 
areas of P n . 

Condition (iii) encourages that shape set S to be minimal ; 
i.e., S will not use any extra shapes to form a complete sum- 
mary of the boundaries of classifier P n . A complete summary 
S (i.e., one satisfying (i) and (ii)) remains a complete summary 
if one of its shapes S G S is added to S either exactly or after 
a small perturbation. In fact, adding a small perturbation S 
of S may actually improve completeness slightly since S can 
be even closer to some high entropy points than S. And if 
S were a good summary, then so too would S. We need the 
minimality condition (iii) to be able to obtain the simplest (i.e., 
smallest) shape set that is a complete summary of the classifier 
boundaries. 

3) Statistical Modeling: The shape set posterior is 

p/ cl V } — p(Y ICApfCA 

P\S \X n ) — p(x ) P\Xn\S)P\S)- 

We build the posterior model P(S\X n ) by modeling the 
likelihood P(X n \S) and the shape set prior P(S). In the 
posterior P(S\X n ) oc P(X n \S)P(S), we will model the 
likelihood P(X n \S) to encourage completeness and the prior 
P(S) to encourage distance between shapes and therefore 


minimality. It makes sense that the data likelihood accounts for 
completeness because completeness requires observed points 
to be close to a shape and the observed points arise from 
the ground truth shapes with the addition of noise. We will 
encourage good summary using a Bayesian loss function that 
grows with increasing distance of the shapes to the point set. 
Finally, we determine the number of shapes l by minimizing 
the expected posterior loss. 

a) Likelihood: Our likelihood will encourage complete- 
ness. For the completeness condition (ii), we are interested 
in making the average squared distance D Xn s °f boundary 
points in X n = {xi , . . . , x n } to shapes in S small: 

^2 _E x eX n d X n ,s( X ) _T,j = l d X n ,s( X i) 

\x n \ \x n \ ’ ( } 

where 

d x n ,s( x ) = ™n\\x - s\\l ( 2 ) 

is the minimum squared distance of a high entropy point x to 
a point on any shape in the collection S = (Si, . . . , Si). 

An observed point x 3 G X n is assumed to have been 
generated from a shape S Zj , where z 3 gives the shape number 
that explains Xj. Given Zj, we model the likelihood of Xj as a 
decreasing function of the minimum distance from Xj to S Zj . 
The closer Xj is to shape S Zj9 the higher the likelihood of 
Xj. The observations Xj are assumed to be independent and 
modeled as 

x j = s i + £ 3 = S j + r 3 n v r j ~ m <7r)> 

where nj is a unit normal to S Zj at Sj and r j = (xj — Sj) • nj . 
Here the noise vector Sj = r 3 n 3 is along a unit normal nj to 
the shape S Zj at the closest shape point Sj to Xj. The scalar 
residual rj is the signed distance along nj from the shape S Zj 
to Xj. We model the observation error Sj by modeling the 
signed residual as a ^(O,^ 2 ) random variable. 

Note that the squared residual r 2 is just the minimum 
distance squared from Xj to the closest point Sj on shape S Zj : 

r 2 , = min \\xj — silo, 

J ses Zj 11 3 M2J 

where the minimum occurs at s — Sj. Let Z = (zi , . . . , z n ). 
Assuming independence of points and that Xj depends only 
on shape S Zj , then P(X n \Z,S ) = Yl] =1 P(xj\zj,S Zj ) = 
YYj=i N{rj |0, cr 2 ). Since rj ~ 7V(0, cr 2 ), it follows that 

P(X n \Z, S) Kcjf n exp E \\ x 3 - s j\\^j 5 ( 3 ) 

for a constant K. Note that if the observed point set X n is close 
to the shapes in S , then P{X n \Z,S) is high. This statement 
assumes, of course, that the correct shape S Zj explaining each 
point Xj has also been identified. 

We can obtain the likelihood P(X n \S) by model- 
ing Z\S and integrating out Z as in P(X n \S) = 
frP{X n \Z,S)P(Z\S)dZ. We could, for example, model 
Z]S by modeling a count vector C = (ci,...,q) which 
holds the number of observations q explained by shape Si. 
Here c t = Y^=i 1 zj=i • We can encourage good summary 
by modeling C ~ multinomial(n, (!//,!//,. ..,!//)) where 



each of the Z shapes in S has the same probability l/l of 
generating an observed point. This would make shape sets with 
any shapes that are from the data quite unlikely because we 
would expect to see points around each shape according to the 
given multinomial distribution. 

It is difficult, however, to optimize over shape sets with 
the hidden random variables Z in our models. Instead, we 
make a simple but accurate and effective approximation in 
our models and assume that the shape S Zj that explains 
observation Xj is the shape in S which is closest to Xj. Thus 
we replace the minimization in equation (3) over Sj G S Zj 
with a minimization sj G S over the entire shape set to obtain 
the approximation 

P(X n \S) = Ka~ n exp ^-^2 I \ x i ~ s j\\^j • ( 4 ) 

From equations (1),(2), we can see that the inner sum in 
equation (4) is just a scaled version \X n \D XnS of our 

completeness measure. We can easily write our likelihood in 

2 

terms of the completeness measure D Xn s . To do so cleanly, 
define a c 2 omplete = a 2 /\X n \. Then 

P(X n \S) = A-a^ iplcte exp (-—fi D 2 XntS ) , 

\ ^complete J 

where another constant factor has been absorbed into K. 


b) Shape Set Prior: We build the shape set prior P(S) 
based on the distances of points on each shape Si to the rest 
of the shape set S-i = S\{Si}. To keep shapes apart from 
one another, we want a large average squared distance from 
points on each shape to the rest of the shapes. Let d 2 s . s ( s i) 
be the minimum squared distance of a point Si G Si to another 
shape Sj: 

(•''•/) = min b ' Sj\\l 

s j ^^3 

Then the squared distance of Si G Si to the shape set S-i is 

d s u S-M) = c mm ds itSj (8i), 


which finds the closest point in the rest of the shapes S-i to 
Si G Si. Finally we average the inter-shape squared distances 
over all points on all shapes to get 
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To keep the shapes apart a priori, we want D s to be large, 
indicating that on average the inter- shape distance is large. 
Equivalently, 1 / Ds should be small. Therefore we model the 
prior for S using the normal distribution 


S ~ N(D S ; 0, crg hapesim ). 


c) Bayesian Loss: Next we define a Bayesian loss 
function that encourages good summary. We can think of the 
summary condition (i) as requiring a small distance from each 
shape S G S to the set of classifier boundary points X n . Let 
d 2 s x ( s ) denote the squared distance from a shape point s G S 
to the point set X n \ 

d s,x n ( s ) = min \\s - x\\ 2 2 . 

xEX n 


— 2 

We capture the average squared distance D s Xn from the shape 
set S to the input points X n by averaging over all points on 
all shapes in S = (Si , . . . ,Si): 

jp _ D fl =i ^2ses a ds a ,x n ( s ) 

S,XU ~ El = 1 \Sa\ 

We define our Bayesian loss function as 

2 

loss (5, 2f n ) = A summar y_D£ , Xn 

The smaller the distance from each shape in S to the point 
set X n , the smaller the loss. Thus minimizing the loss will 
encourage good summary. 

4) Shape Fitting Method: Our shape fitting method has two 
main steps: 

Step 1 Minimize the expected posterior loss 

g(l) = E[\oss(S,X n )}, \S\ = l 

over l to obtain the number of shapes Z* 

Step 2 Compute the MAP shape set S* ,L for sets of size 

P 

a) Determining the Number of Shapes: We assume that 
we can apriori limit the number of shapes l to some set C. For 
example, if we know that there will not be more than five 
boundaries then we can set C = {1, 2, 3, 4, 5}. 

For each Z G £, we compute the expected posterior loss 

g(l) = E[\oss(S,X)} = [ \oss(S,X n )P(S\X n )dS. 

Here we denote the shape set posterior probability distribution 
for shape sets with a fixed number of shapes as P(S\X n ). Then 
we choose the number of shapes to minimize the expected 
posterior loss: 

Z* = argming(Z). 

5) Our Shape Set Posterior Sampling Method: For a fixed 
shape set size \S\ = Z, we will draw samples from the posterior 
P(S \X n ) oc P(X n \S)P(S) using an iterative procedure. 
Shape set samples S with a small value for 

mog(P(X n \S)P(S)) = -log(P(X n \S)) - log (P(S)) 

should be more likely to occur. 

V. Experiments and Results 
A. Evaluation on Artificial Data Sets 

1) Active Learning: We illustrate the behavior of our 
approach using a 2D artificial data set and a quadratic boundary 
function normalized to a unit square. Starting with a low 
number of N init = 126 randomly selected initial data points, 
the active learning procedure selects N = 500 new data points 
according to different candidate selection rules (random, ALC, 
ALM, El, and boundary El). N has been selected this large 
for visualization purposes. Figure 3 shows, how the different 
selection algorithms behave. Our goal is to find many data 
points near the threshold curve in order to enable accurate 
representation and to facilitate subsequent shape estimation. 




Fig. 3: Candidate points during active learning: (A) random 
selection, (B) ALC, (C) ALM, (D) El, and (E) boundary-EI. 
Circles: initial data points. Solid: points added during active 
learning (colored according to experiment outcome). 


near the threshold (high closeness) but extremely low coverage. 
Figure 4 shows this trade space. For each of the newly added 
points, we calculate d as the minimal distance of that point to 
the boundary. Obviously, small values should be preferred, as 
such points close to the threshold help to accurately estimate its 
shape. Figure 4 A shows a histogram of distances d for various 
update rules. Whereas random and ALC have large numbers 
of points that are far away from the threshold surface, ALM 
seems to perform best for this metric. However, Figure 4B 
reveals that ALM only covers a very small portion of the 
threshold surface. Random selection provides the best coverage 
here. With our analysis goal in mind, our boundary- aware El 
metric features a good overall coverage and a high density of 
points close to actual threshold surface. 

Our boundary metric is parameterized by the parameter a 
(see Section III-C). This parameter influences the width of the 
’’band” around the threshold surface that is considered for the 
selection of the candidate point. Figure 5 shows runs with 
several values of a. It seems that values around a = 0.8 
produce the best results; values of a that are too small or 
too large tend to lead to a situation, where the new points are 
located too far from the threshold surface. 



Fig. 5: Boundary-EI for a = 0.2, 0.5, 0.8, 1 


On the other hand, the entire area should be considered as 
well in order not to miss any other boundary. 

The random Monte-Carlo style selection (Figure 3A) needs 
a prohibitively large N for reasonable results. The classical 
ALC [13] (Figure 3B) finds many points near the boundary, 
but still too many data points are away from the curve, 
demanding large N. Other algorithms are too localized and 
do not even explore the entire threshold curve (Figure 3C, 
D). Our approach (Figure 3E) tries to find a suitable balance 
between both requirements. 



Fig. 4: A: Histogram for distances d of candidate points 
from boundary. Leftmost bars cropped for better visibility. B: 
Histograms of boundary coverage. 


There is a trade-space between closeness of selected points 
to the safety boundary and a good curve coverage. For exam- 
ple, a greedy algorithm might always select the same point 


The performance of the active learning method and the 
shape estimation can also be assessed by analyzing the al- 
gorithm convergence, i.e., how many new data points are 
necessary, before the estimated shape parameters are close to 
the ground truth that has been used to generate the artificial 
data set. For e xample, for a single hyperplane boundary, we 

obtain C = \J \9\ — \9\, where 6 and 6 are the ground truth 
and the estimated parameters, respectively. Figure 6 compares 
the convergence of random selection, ALC, and our method 
over 10 runs and shows a superior performance of our method. 





Fig. 6: Convergence C for random (red), ALC (black), and 
boundary-EI (green) over learning iterations. 

2) Shape Selection: To assess the performance of our shape 
selection and estimation method, we carried out experiments 
with boundaries in the shape of hyperspheres. Table I A shows 
the results in the two-dimensional case. With a total 126+700 
data points, sphere Si was correctly recognized in all 25 
runs; S 2 was only correctly estimated in 5 of the runs. The 


table shows the ground-truth values, means and variances for 
the successful runs. 74 data points were selected for shape 
estimation. 


TABLE I: Parameters for 2D (A) and 5D spheres (B). 



true value 

AO ) 

Xo 

0.3 

0.295(1.4e-5) 

yo 

0.3 

0.289(3e-5) 

ro 

0.2 

0.20(6.6e-6) 

X 1 

0.7 

0.715(8.5e-5) 

y i 

0.7 

0.72(8.6e-5) 

ro 

0.2 

0.20(5.2e-5) 



true value 

AO ) 

Cl 

0.3 

0.29(7e-3) 

C2 

0.3 

0.26(5e-3) 

C3 

0.3 

0.32(8e-3) 

C4 

0.3 

0.31(7e-3) 

C5 

0.3 

0.27(9e-3) 

r 

0.3 

0.29(8e-4) 


Table IB shows the situation in a 5D space. With the 
centers located at c\ = (0.3, 0.3, 0.3, 0.3, 0.3) T , and c\ = 
(0.7, 0.7, 0.7, 0.7, 0.7) t , respectively and radius r = 0.3. 
Active learning selected 1000 data points, 155 of which were 
selected for shape estimation. Here, the results are much worse. 
E.g., the second sphere was not recognized in any of the 10 
runs, indicating that for a 5D space, the number of data points 
must be considerable larger. Future work will address this 
issue. 


Fig. 7: A: Actual boundary (blue) and estimated boundary 
(green) over w p ,w q , Ki at . B: Estimated hypersphere shape. 

TTSAFE are the aircraft trajectories as well as a number 
of configuration parameters. During each run, the number of 
detected losses of separation (“numlos”) as well as the time to 
loss of separation (“ttlos”) is reported. For a given scenario, 
low numbers of numlos and late detection (small ttlos) mean 
unsafe situations. We analyzed the behavior of the TTSAFE 
system with respect to different parameter values and altitude 
biases. Here we focus on the altitude biases, which can be 
caused by noisy or faulty radar measurements. 




B. IFCS Data Set 

The Intelligent Flight Control System (IFCS) is a damage- 
adaptive Neural Networks (NN) based flight control system 
developed by NASA and test-flown on a manned F-15 aircraft 
[2]. An on-line trained NN provides control augmentation to 
dynamically counteract damages to the aircraft. For our exper- 
iments, we considered this system as a black box, controlled 
by numerous parameters (e.g., NN weights, controller gains, or 
learning rate). A simulation run was considered to be success- 
ful, if, after an injected damage, the aircraft remained stable 
for at least 20 seconds. After an initial parameter sensitivity 
analysis, we selected the parameters w p ,w q ,w r , Ki at , and ( 
for further analysis, where the Wi are proportional gains of 
the controllers, Ki at the lateral stick gain, and ( a damping 
coefficient. We generated a combinatorial data set of 32,768 
data points, out of which 7,992 runs were successful. 

A boundary over these parameters exist in a shape of a 
hypersphere. This spherical shape is a consequence of the 
IFCS design, and the shape can be described by ( Wp ^° ) 2 + 

r«A 2 ~ y ° ) 2 = r--»»~* 0 ) 2 = C x fa - K^. This stability 
boundary is is parameterized by unknown fa . , yo , zq are 

design-time constants. 

Figure 7 A shows the actual and estimated boundary in a 
projection into w p ,w q , and Ki at . For our shape fitting and 
estimating experiment, we used 1000+5000 data points. The 
shape parameters for the boundary in Figure 7B were estimated 
based upon 485 points near the boundary within an e-band of 
width 0.2. 



Fig. 8: A: Active learning of boundaries. B: Different values 
of numlos over different altitude biases for aircraft 1 and 2. 

Figure 8 A shows a boundary with respect to the TTSAFE 
parameters “MinHorizSepNmi” (minimal horizontal separation 
in Nautical miles) and how often the detection algorithm is 
executed “CheckPeriodMinutes” in a normalized space. Red 
and green circles denote the initial data points; triangles the 
points added by active learning. The blue line visualizes the 
estimated (hyperplane) boundary. A bias in the altitudes of 
the two converging aircraft can confuse TTSAFE: the aircraft 
might appear closer to each other than they actually are, 
causing a false alarm. Alternatively, they might appear farther 
apart than in reality. This is a dangerous situation as TTSAFE 
would not detect the loss of separation. In Figure 8B we show 
the space spanned by the two bias parameters. Blue dots are 
in a safe region (higher numbers of numlos); red and magenta 
dots represent unsafe regions. 

Figure 9A shows results of the above experiment with 
respect to time to loss (ttlos). It turned out that this boundary 
is rather small and hidden in the entire parameter space 
(Figure 9B). 


C. Experiments with TTSAFE 

The Terminal Tactical Separation Assured Flight Envi- 
ronment (TTSAFE) is a software tool to predict violations 
of separation rules between aircraft near an airport [17]. 
TTSAFE uses novel algorithms for predicting the trajectories 
and encodes all complex terminal separation rules. Thus, we 
decided to analyze TTSAFE as a black box. The inputs to 


VI. Conclusion and Discussion 

In this paper, we described a statistical modeling ap- 
proach for the online detection and characterization of safety 
boundaries. Given a library of candidate shapes, location and 
parameters of the boundaries are estimated using an advanced 
Bayesian modeling approach. The results of our boundary 
analysis is provided in a form understandable by the domain 



[15] P. Ranjan, D. Bingham, and G. Michailidis, “Sequential experiment 
design for contour estimation from complex computer codes,” Techno- 
metrics, vol. 50, no. 4, pp. 527-541, 2008. 

[16] D. Jones, M. Schonlau, and W. J. Welch, “Efficient global optimization 
of expensive black box functions,” Journal of Global Optimization, 
vol. 13, pp. 455-492, 1998. 

[17] H. Tang, J. Robinson, and D. Denery, “Tactical conflict detection in 
terminal airspace,” in 10th AIAA Aviation Technology, Integration, and 
Operations (ATIO) Conference, 2010. 


Fig. 9: A: Boundary of ttlos over altitude biases. New data 
points added by active learning with ttlos > 0.5 (blue) and 
ttlos <0.5 (red). B: Boundary shown in A (red) in otherwise 
insensitive parameter space. 


expert. Our active learning procedure uses a boundary-aware 
metric to quickly and effectively find new data points near the 
boundary. Experiments show that existing boundaries could 
be reliably found and their shape parameters estimated with 
confidence interval. Our boundary-EI substantially improved 
the selection of new data points. 

In our future work, we will investigate whether a tighter 
coupling of the selection metric and the posteriors of shape 
detection can improve the active learning performance. 
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